function f = bifConstant(p0, gammaLower, gammaUpper)

g = 9.8;
m = 10000;
h = (gammaUpper-gammaLower)/m;
for i = 1:m
	gamma(i) = (i-1)*h+ gammaLower;
	f(i) = tanh(sqrt(2*abs(p0)/gamma(i))) - ...
		2*abs(p0)*gamma(i)/(g +gamma(i)*sqrt(2*abs(p0)*gamma(i)));
end
plot(gamma, f);
title({['p_0 = ', num2str(p0)]});
end


